Built with R version 3.6.3

0.1 Setup

0.1.1 Environment

knitr::opts_chunk$set(
  echo = TRUE,
  message = FALSE,
  warning = FALSE,
  eval = TRUE, 
  tidy = TRUE
)

# Environment Set up
rm(list = ls()) #Clean workspace
cat("\014")     #Clean Console

gc() # Free memory

##          used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
## Ncells 487947 26.1    1073392 57.4         NA   641482 34.3
## Vcells 926847  7.1    8388608 64.0      32768  1768755 13.5

###################
# Install packages
###################

pkgs <- c("remedy", "Seurat", "dplyr", 
          "rstudioapi", "cowplot", "ggplot2", 
          "grid", "gridExtra", "styler", "stringr", "msigdbr")


for(i in 1:length(pkgs)){
  if(!require(pkgs[i], character.only = T)){
    install.packages(pkgs[i])
    require(pkgs[i], character.only = T)
  }else{
    require(pkgs[i], character.only = T)
  }
}

pkgs <- c("gplots", "fgsea", "biomaRt",
          "clusterProfiler", "GSEABase", "Bitr",
          "org.Mm.eg.db")

for(i in 1:length(pkgs)){
  if(!require(pkgs[i], character.only = T)){
    BiocManager::install(pkgs[i])
    require(pkgs[i], character.only = T)
  }else{
    require(pkgs[i], character.only = T)
  }
}

## Warning in library(package, lib.loc = lib.loc, character.only = TRUE,
## logical.return = TRUE, : there is no package called 'Bitr'

## Warning: package 'Bitr' is not available (for R version 3.6.3)

## Warning in library(package, lib.loc = lib.loc, character.only = TRUE,
## logical.return = TRUE, : there is no package called 'Bitr'

####################
# Github packages
####################

# "EnhancedVolcano" package from Github repo
# devtools::install_github('kevinblighe/EnhancedVolcano')
library(EnhancedVolcano)

# Colour scheme
heatmap.blue.col <- "#2B4FA2"
heatmap.red.col <- "#ED1C24"
white.col <- "#FFFFFF"

DNAM1.neg.col <- "#FF0000" # Red 
DNAM1.dim.col <- "#A0A0A4" # Grey
DNAM1.high.col <- "#0000C0" # Blue
col.scheme <- c(DNAM1.neg.col, DNAM1.dim.col, DNAM1.high.col)

# Set working directory to source file location
#setwd(dirname(getActiveDocumentContext()$path))
#working.dir <- getwd()


# create output directories
if(!dir.exists("output")){dir.create("output", recursive = T)}
if(!dir.exists("output/figures")){dir.create("output/figures", recursive = T)}
if(!dir.exists("output/tables")){dir.create("output/tables", recursive = T)}
if(!dir.exists("output/QC")){dir.create("output/QC", recursive = T)}

0.1.2 Reading data

# Load dataset
DNAM1.data <- Read10X(data.dir = "Input_data/10X/")

# Initialize the Seurat object with the raw (non-normalized data).  Keep all
# genes expressed in >= 3 cells.  Keep all cells with at least 200 detected genes
seurat.object <- CreateSeuratObject(counts = DNAM1.data$`Gene Expression`, min.cells = 3, 
    min.features = 200)

seurat.object  # 13,172 genes and 1,907 cells

## An object of class Seurat 
## 13172 features across 1907 samples within 1 assay 
## Active assay: RNA (13172 features, 0 variable features)

0.2 QC and Data filtering

0.2.1 QC & Normalisation

# Get mitochondria genes and calculate percentage of genes
mito.genes <- grep(pattern = "^mt-", x = rownames(x = seurat.object@assays$RNA@data), 
    value = TRUE)
mito.genes

##  [1] "mt-Nd1"  "mt-Nd2"  "mt-Co1"  "mt-Co2"  "mt-Atp8" "mt-Atp6" "mt-Co3" 
##  [8] "mt-Nd3"  "mt-Nd4l" "mt-Nd4"  "mt-Nd5"  "mt-Nd6"  "mt-Cytb"

percent.mito <- Matrix::colSums(seurat.object@assays$RNA@data[mito.genes, ])/Matrix::colSums(seurat.object@assays$RNA@data)

seurat.object <- AddMetaData(object = seurat.object, metadata = percent.mito, col.name = "percent.mito")


# Visualise mitochondria percentage
VlnPlot(object = seurat.object, features = c("nFeature_RNA", "nCount_RNA", "percent.mito"))

dev.copy(pdf, "output/QC/QC_nFeature_nCount_RNA_percentMito.pdf")

## pdf 
##   3

dev.off()

## quartz_off_screen 
##                 2

head(seurat.object@meta.data)

str(seurat.object@meta.data)

## 'data.frame':    1907 obs. of  4 variables:
##  $ orig.ident  : Factor w/ 1 level "SeuratProject": 1 1 1 1 1 1 1 1 1 1 ...
##  $ nCount_RNA  : num  9729 8060 634 7692 2296 ...
##  $ nFeature_RNA: int  2510 2117 411 2230 1045 2140 2589 1564 2767 1849 ...
##  $ percent.mito: num  0.0363 0.0362 0.0994 0.0261 0.2134 ...

################### Filter cells

# Remove cells with < 200 genes, and remove cells with > 0.05% mitochondrial
# genes

seurat.object.filt <- subset(seurat.object, subset = nFeature_RNA > 200 & nFeature_RNA < 
    6000 & percent.mito < 0.05)

seurat.object.filt  #13,172 genes and 1,526 cells remain

## An object of class Seurat 
## 13172 features across 1526 samples within 1 assay 
## Active assay: RNA (13172 features, 0 variable features)

# Vis filtering

# Before filtering
VlnPlot(object = seurat.object, features = c("nFeature_RNA", "nCount_RNA", "percent.mito"))

# After filtering
VlnPlot(object = seurat.object.filt, features = c("nFeature_RNA", "nCount_RNA", "percent.mito"))

dev.copy(pdf, "output/QC/QC_nFeature_nCount_RNA_percentMito_Post_filtering.pdf")

## pdf 
##   3

dev.off()

## quartz_off_screen 
##                 2

seurat.object  # 13,172 genes across 1,907 cells

## An object of class Seurat 
## 13172 features across 1907 samples within 1 assay 
## Active assay: RNA (13172 features, 0 variable features)

seurat.object.filt  # 13,172 genes across 1,526 cells

## An object of class Seurat 
## 13172 features across 1526 samples within 1 assay 
## Active assay: RNA (13172 features, 0 variable features)

# Therefore 381 cells were removed


# Normalise RNA data by logNormalize
seurat.object.filt <- NormalizeData(object = seurat.object.filt, normalization.method = "LogNormalize", 
    scale.factor = 10000, verbose = TRUE)


# Find variable genes
seurat.object.filt <- FindVariableFeatures(object = seurat.object.filt, verbose = TRUE)

# Scale and regress variables
seurat.object.filt <- ScaleData(object = seurat.object.filt, features = rownames(seurat.object.filt), 
    vars.to.regress = c("nCount_RNA", "percent.mito"), model.use = "linear", block.size = 2000, 
    do.scale = TRUE, do.center = TRUE, verbose = TRUE)

0.2.2 Dimension reduction

# Run PCA
seurat.object.filt <- RunPCA(seurat.object.filt, verbose = FALSE)

ElbowPlot(seurat.object.filt, ndims = 50)

dev.copy(pdf, "output/figures/PC_elbowplot.pdf")

## pdf 
##   3

dev.off()

## quartz_off_screen 
##                 2

# elbow @ ~20

PCAPlot(seurat.object.filt, dims = c(1, 2)) + NoLegend()

dev.copy(pdf, "output/figures/PCA_dims1_2.pdf")

## pdf 
##   3

dev.off()

## quartz_off_screen 
##                 2

seurat.object.filt <- FindNeighbors(seurat.object.filt, dims = 1:20)

seurat.object.filt <- FindClusters(seurat.object.filt, resolution = 0.1)

## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 1526
## Number of edges: 55479
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9179
## Number of communities: 3
## Elapsed time: 0 seconds

########## tSNE

# tSNE perplex picked based on looped running of RunTSNE at differing values

seurat.object.filt <- RunTSNE(seurat.object.filt, reduction = "pca", dims = 1:20, 
    perplexity = 30, method = "FIt-SNE", seed.use = 42, dim.embed = 2)

DimPlot(seurat.object.filt, reduction = "tsne") + ggtitle(paste0("Perplexity set at ", 
    30))

dev.copy(pdf, "output/figures/tSNE_all_cells_res01_clusts.pdf")

## pdf 
##   3

dev.off()

## quartz_off_screen 
##                 2

######### UMAP

# UMAP min.dist and n.neighbors values were picked based on looped running of
# RunUMAP across a range of values

seurat.object.filt <- RunUMAP(object = seurat.object.filt, reduction = "pca", dims = 1:20, 
    umap.method = "uwot", n.neighbors = 20, min.dist = 0.5, seed.use = 42)

DimPlot(object = seurat.object.filt, reduction = "umap", label = TRUE) + ggtitle(paste0("UMAP Min dist = ", 
    0.5)) + theme(plot.title = element_text(hjust = 0.5))

dev.copy(pdf, "output/figures/UMAP_all_cells_res01_clusts.pdf")

## pdf 
##   3

dev.off()

## quartz_off_screen 
##                 2

0.2.3 Adding ADT data and normalising/filtering

################################## ADD ADT data to seurat object Data needs to be transposed Additionally, cells
################################## filtered during mRNA QC need to be removed

ADT.data <- as.matrix(DNAM1.data$`Antibody Capture`)

dim(ADT.data)

## [1]    1 1936

length(rownames(seurat.object.filt@meta.data))

## [1] 1526

# Determine which cells remain in seurat dataset and keep only these from ADT
# dataset
cell.ids.keep <- rownames(seurat.object.filt@meta.data)
keep.logic <- colnames(ADT.data) %in% cell.ids.keep
ADT.data <- ADT.data[, keep.logic]
ADT.data <- t(as.matrix(ADT.data))
rownames(ADT.data) <- "DNAM1"

# Add ADT data to seurat object
seurat.object.filt[["ADT"]] <- CreateAssayObject(counts = ADT.data)
seurat.object.filt

## An object of class Seurat 
## 13173 features across 1526 samples within 2 assays 
## Active assay: RNA (13172 features, 2000 variable features)
##  1 other assay present: ADT
##  3 dimensional reductions calculated: pca, tsne, umap

# ADT Normalisation
seurat.object.filt <- NormalizeData(seurat.object.filt, assay = "ADT", normalization.method = "CLR")
# ADT scaling
seurat.object.filt <- ScaleData(seurat.object.filt, assay = "ADT")

0.2.4 Subset data by ADT_Cd226 expression

############################################### Subset cells based on Cd226/DNAM1 ADT expression

# Cells are subsetted by ADT (Cd226/DNAM1) expression into Neg, Dim, High
# populations DNAM1 Neg = ~25% of cells with the lowest expression DNAM1 Dim =
# next ~50% of cells DNAM1 High = top ~25% of cells with highest expression


seurat.object.filt  # 1,526 samples/cells

## An object of class Seurat 
## 13173 features across 1526 samples within 2 assays 
## Active assay: RNA (13172 features, 2000 variable features)
##  1 other assay present: ADT
##  3 dimensional reductions calculated: pca, tsne, umap

# Get quantile cutoff values
q.vals <- quantile(t(as.matrix(seurat.object.filt@assays$ADT@data)), probs = c(0.25, 
    0.75))
print(q.vals)

##       25%       75% 
## 0.3774276 1.2486700

# Check number of cells in each subset
sum(seurat.object.filt@assays$ADT@data <= q.vals[1])  # 384

## [1] 384

sum(seurat.object.filt@assays$ADT@data > q.vals[1] & seurat.object.filt@assays$ADT@data < 
    q.vals[2])  # 756

## [1] 756

sum(seurat.object.filt@assays$ADT@data >= q.vals[2])  # 386

## [1] 386

# Get logical variables of which cells fit in each subset
logic.var.neg <- seurat.object.filt@assays$ADT@data <= q.vals[1]
logic.var.dim <- seurat.object.filt@assays$ADT@data > q.vals[1] & seurat.object.filt@assays$ADT@data < 
    q.vals[2]
logic.var.high <- seurat.object.filt@assays$ADT@data >= q.vals[2]

# Get cell IDs for each subset
DNAM.neg <- colnames(seurat.object.filt@assays$ADT@data)[logic.var.neg]
DNAM.dim <- colnames(seurat.object.filt@assays$ADT@data)[logic.var.dim]
DNAM.high <- colnames(seurat.object.filt@assays$ADT@data)[logic.var.high]

# Create a vector with new subset annotations
cell.ids <- rownames(seurat.object.filt@meta.data)
neg.logic <- cell.ids %in% DNAM.neg
dim.logic <- cell.ids %in% DNAM.dim
high.logic <- cell.ids %in% DNAM.high

DNAM.clust <- cell.ids
DNAM.clust[neg.logic] <- "Neg"
DNAM.clust[dim.logic] <- "Dim"
DNAM.clust[high.logic] <- "High"
DNAM.clust

##    [1] "Neg"  "Dim"  "High" "High" "Neg"  "Neg"  "Dim"  "Dim"  "Neg"  "Neg" 
##   [11] "Dim"  "Dim"  "Dim"  "High" "Dim"  "Dim"  "High" "Neg"  "Dim"  "Dim" 
##   [21] "Neg"  "Dim"  "Dim"  "Neg"  "Neg"  "High" "High" "Neg"  "Dim"  "High"
##   [31] "High" "Dim"  "Neg"  "High" "High" "Dim"  "High" "Neg"  "Dim"  "Neg" 
##   [41] "Dim"  "High" "Neg"  "Neg"  "Dim"  "High" "Neg"  "Dim"  "Dim"  "Dim" 
##   [51] "Dim"  "Dim"  "Neg"  "Dim"  "Dim"  "Dim"  "Dim"  "Neg"  "Dim"  "Dim" 
##   [61] "Dim"  "Dim"  "Dim"  "Dim"  "Dim"  "High" "Dim"  "Neg"  "Neg"  "Neg" 
##   [71] "Neg"  "Dim"  "Dim"  "High" "Dim"  "Dim"  "Neg"  "High" "High" "High"
##   [81] "High" "High" "High" "Dim"  "Dim"  "Dim"  "Neg"  "High" "High" "Dim" 
##   [91] "Dim"  "Dim"  "Dim"  "Dim"  "Dim"  "Dim"  "Neg"  "Dim"  "High" "High"
##  [101] "High" "Dim"  "Dim"  "High" "Dim"  "Dim"  "High" "Dim"  "Dim"  "Dim" 
##  [111] "Dim"  "High" "Dim"  "Neg"  "Neg"  "Dim"  "Dim"  "High" "Dim"  "Dim" 
##  [121] "High" "Dim"  "High" "High" "Dim"  "Dim"  "Dim"  "High" "Neg"  "High"
##  [131] "Neg"  "Dim"  "Dim"  "Dim"  "Dim"  "Dim"  "Neg"  "Dim"  "High" "Dim" 
##  [141] "Dim"  "Dim"  "High" "Neg"  "High" "Dim"  "Dim"  "Neg"  "Neg"  "Neg" 
##  [151] "Dim"  "Dim"  "High" "Neg"  "High" "Neg"  "Dim"  "Dim"  "Neg"  "Neg" 
##  [161] "High" "Dim"  "High" "Neg"  "Dim"  "Neg"  "Dim"  "Neg"  "Dim"  "Dim" 
##  [171] "Dim"  "Dim"  "Dim"  "Dim"  "Dim"  "Neg"  "High" "Dim"  "High" "Neg" 
##  [181] "High" "Neg"  "Dim"  "Dim"  "Neg"  "High" "High" "Neg"  "High" "Dim" 
##  [191] "Neg"  "Dim"  "Dim"  "High" "Neg"  "High" "Neg"  "Neg"  "Neg"  "Neg" 
##  [201] "High" "Dim"  "Dim"  "Neg"  "High" "Dim"  "Dim"  "Dim"  "Dim"  "High"
##  [211] "Neg"  "Neg"  "Dim"  "Dim"  "Dim"  "Dim"  "Dim"  "Neg"  "Dim"  "Dim" 
##  [221] "Neg"  "High" "Dim"  "Dim"  "Neg"  "Dim"  "High" "High" "Neg"  "High"
##  [231] "Dim"  "High" "Dim"  "High" "Dim"  "Neg"  "High" "Dim"  "High" "Dim" 
##  [241] "Neg"  "Neg"  "Dim"  "Dim"  "Dim"  "Dim"  "High" "Dim"  "Dim"  "Dim" 
##  [251] "High" "Dim"  "Dim"  "High" "Dim"  "Dim"  "Neg"  "Neg"  "Neg"  "Dim" 
##  [261] "High" "Neg"  "Dim"  "Dim"  "Dim"  "Dim"  "Neg"  "High" "Neg"  "High"
##  [271] "Dim"  "Neg"  "Neg"  "High" "Neg"  "Neg"  "Dim"  "High" "High" "Neg" 
##  [281] "Neg"  "Dim"  "High" "Neg"  "High" "Dim"  "High" "High" "Dim"  "Dim" 
##  [291] "Dim"  "High" "High" "Dim"  "Dim"  "High" "Dim"  "High" "High" "Dim" 
##  [301] "Neg"  "Dim"  "Dim"  "Dim"  "High" "Dim"  "Neg"  "Neg"  "Neg"  "High"
##  [311] "High" "Dim"  "Neg"  "Neg"  "Dim"  "Dim"  "Dim"  "High" "Dim"  "Dim" 
##  [321] "Dim"  "Dim"  "Neg"  "Neg"  "Dim"  "Dim"  "High" "Dim"  "Neg"  "Dim" 
##  [331] "High" "Neg"  "Dim"  "High" "High" "High" "High" "Dim"  "Dim"  "Neg" 
##  [341] "Dim"  "High" "Dim"  "Dim"  "High" "High" "Dim"  "Neg"  "Dim"  "Dim" 
##  [351] "Dim"  "Dim"  "Dim"  "Dim"  "Dim"  "Dim"  "Dim"  "Dim"  "High" "Neg" 
##  [361] "High" "Neg"  "Dim"  "Dim"  "Neg"  "Neg"  "Dim"  "Dim"  "Dim"  "Dim" 
##  [371] "Dim"  "Neg"  "Dim"  "High" "Dim"  "Dim"  "Dim"  "Dim"  "Dim"  "Neg" 
##  [381] "High" "Dim"  "Dim"  "Dim"  "High" "Dim"  "High" "Neg"  "Neg"  "Dim" 
##  [391] "Neg"  "Dim"  "Dim"  "High" "Dim"  "Neg"  "Neg"  "Dim"  "High" "Neg" 
##  [401] "Neg"  "Neg"  "Neg"  "Dim"  "Dim"  "Neg"  "Neg"  "High" "Dim"  "Dim" 
##  [411] "Dim"  "Dim"  "Neg"  "Neg"  "Dim"  "High" "Dim"  "High" "High" "High"
##  [421] "Dim"  "High" "Dim"  "Dim"  "Dim"  "Neg"  "Dim"  "High" "Dim"  "Neg" 
##  [431] "High" "Dim"  "Dim"  "Dim"  "Neg"  "Dim"  "Dim"  "Dim"  "Dim"  "High"
##  [441] "Dim"  "High" "High" "Neg"  "High" "Dim"  "High" "Neg"  "Dim"  "Dim" 
##  [451] "Neg"  "Dim"  "Neg"  "Dim"  "Dim"  "Dim"  "High" "Neg"  "High" "High"
##  [461] "Dim"  "Neg"  "High" "High" "Dim"  "Neg"  "Dim"  "Dim"  "Neg"  "Dim" 
##  [471] "Neg"  "Dim"  "Neg"  "High" "Dim"  "Neg"  "Dim"  "Dim"  "Dim"  "Dim" 
##  [481] "Dim"  "Dim"  "Dim"  "Dim"  "High" "High" "Neg"  "High" "Dim"  "Dim" 
##  [491] "Dim"  "Dim"  "Neg"  "Neg"  "High" "Dim"  "High" "Neg"  "Dim"  "Dim" 
##  [501] "Dim"  "Neg"  "High" "High" "Dim"  "Dim"  "High" "Dim"  "Neg"  "Neg" 
##  [511] "High" "Dim"  "Dim"  "High" "High" "Dim"  "Dim"  "Dim"  "Dim"  "Dim" 
##  [521] "High" "Neg"  "Dim"  "Neg"  "Neg"  "Dim"  "Neg"  "Dim"  "Dim"  "High"
##  [531] "Neg"  "Neg"  "Dim"  "Dim"  "Dim"  "Dim"  "High" "Dim"  "Neg"  "High"
##  [541] "Dim"  "Dim"  "High" "High" "Dim"  "Dim"  "High" "High" "Dim"  "Dim" 
##  [551] "Dim"  "High" "Dim"  "Dim"  "Dim"  "High" "Dim"  "High" "Dim"  "High"
##  [561] "Neg"  "Dim"  "Neg"  "Dim"  "Dim"  "Dim"  "Dim"  "High" "High" "Dim" 
##  [571] "Dim"  "Dim"  "High" "High" "Neg"  "High" "Neg"  "Dim"  "Neg"  "Dim" 
##  [581] "High" "High" "Dim"  "Neg"  "Dim"  "Dim"  "High" "High" "Neg"  "Neg" 
##  [591] "Dim"  "Dim"  "Dim"  "Neg"  "Dim"  "Dim"  "Dim"  "Dim"  "Neg"  "Neg" 
##  [601] "Neg"  "Neg"  "Dim"  "Dim"  "Dim"  "High" "Dim"  "Neg"  "Dim"  "Dim" 
##  [611] "High" "Neg"  "Dim"  "High" "Neg"  "Neg"  "Neg"  "Neg"  "Neg"  "High"
##  [621] "High" "High" "Dim"  "Neg"  "Neg"  "Dim"  "Dim"  "Dim"  "Dim"  "Dim" 
##  [631] "Dim"  "Dim"  "Neg"  "Dim"  "Dim"  "Dim"  "Dim"  "Dim"  "Dim"  "Neg" 
##  [641] "Neg"  "Neg"  "Dim"  "Dim"  "Neg"  "High" "High" "Dim"  "Neg"  "Neg" 
##  [651] "Dim"  "Neg"  "High" "Neg"  "Dim"  "Dim"  "Dim"  "High" "Dim"  "Neg" 
##  [661] "Dim"  "Dim"  "Dim"  "Dim"  "Dim"  "Dim"  "High" "Dim"  "High" "Dim" 
##  [671] "Neg"  "High" "Neg"  "Dim"  "Dim"  "Dim"  "Neg"  "High" "Neg"  "High"
##  [681] "Neg"  "Neg"  "Dim"  "Neg"  "Dim"  "Dim"  "Dim"  "Dim"  "Dim"  "High"
##  [691] "High" "High" "Dim"  "High" "Dim"  "High" "Dim"  "Neg"  "High" "Dim" 
##  [701] "Dim"  "Dim"  "Neg"  "Dim"  "Dim"  "High" "Neg"  "Dim"  "Dim"  "High"
##  [711] "Neg"  "High" "Neg"  "Dim"  "Dim"  "Dim"  "Dim"  "Dim"  "Dim"  "Dim" 
##  [721] "High" "Dim"  "Dim"  "Dim"  "High" "Dim"  "High" "Neg"  "Dim"  "High"
##  [731] "Dim"  "Neg"  "Neg"  "Dim"  "Dim"  "Neg"  "Dim"  "High" "Dim"  "Neg" 
##  [741] "Dim"  "Dim"  "Neg"  "Dim"  "High" "Dim"  "High" "Neg"  "Dim"  "High"
##  [751] "Dim"  "Dim"  "Dim"  "Neg"  "Dim"  "Dim"  "Neg"  "High" "Dim"  "Dim" 
##  [761] "Dim"  "Neg"  "Dim"  "High" "Dim"  "Dim"  "Neg"  "Neg"  "Dim"  "Neg" 
##  [771] "High" "Dim"  "Dim"  "Dim"  "High" "High" "Dim"  "Dim"  "Dim"  "Neg" 
##  [781] "High" "Dim"  "Neg"  "Neg"  "Dim"  "Dim"  "Dim"  "Neg"  "Dim"  "High"
##  [791] "Dim"  "High" "Dim"  "High" "High" "Dim"  "High" "Dim"  "Neg"  "Neg" 
##  [801] "Dim"  "Dim"  "Dim"  "Neg"  "Neg"  "Neg"  "Neg"  "High" "Dim"  "Neg" 
##  [811] "Dim"  "High" "Dim"  "Neg"  "Neg"  "Dim"  "Neg"  "Dim"  "Dim"  "Dim" 
##  [821] "High" "Neg"  "High" "Dim"  "Dim"  "Neg"  "High" "Neg"  "Neg"  "Dim" 
##  [831] "High" "Neg"  "High" "Dim"  "Neg"  "Dim"  "High" "High" "Neg"  "Dim" 
##  [841] "Dim"  "Dim"  "Dim"  "Dim"  "Dim"  "Dim"  "Dim"  "Neg"  "Dim"  "Neg" 
##  [851] "Dim"  "High" "Neg"  "Neg"  "Dim"  "Dim"  "High" "Neg"  "High" "High"
##  [861] "Dim"  "Dim"  "Dim"  "Neg"  "Neg"  "Neg"  "Dim"  "High" "High" "Dim" 
##  [871] "Dim"  "Dim"  "Dim"  "Dim"  "Dim"  "Dim"  "High" "High" "Dim"  "High"
##  [881] "Neg"  "Neg"  "High" "Neg"  "Dim"  "Neg"  "Dim"  "High" "High" "High"
##  [891] "Neg"  "High" "Dim"  "Dim"  "High" "Dim"  "Dim"  "Dim"  "Neg"  "High"
##  [901] "High" "Dim"  "Dim"  "Dim"  "Dim"  "Neg"  "High" "Neg"  "High" "Dim" 
##  [911] "Dim"  "High" "Dim"  "High" "Dim"  "Neg"  "Dim"  "Dim"  "High" "High"
##  [921] "High" "High" "High" "High" "High" "High" "High" "High" "Dim"  "Dim" 
##  [931] "Dim"  "Dim"  "Dim"  "Neg"  "Dim"  "Neg"  "Dim"  "High" "Dim"  "Dim" 
##  [941] "Dim"  "Neg"  "Neg"  "Dim"  "Dim"  "Dim"  "High" "Dim"  "High" "Neg" 
##  [951] "High" "Dim"  "High" "Neg"  "High" "Dim"  "Neg"  "Dim"  "Dim"  "Dim" 
##  [961] "Dim"  "High" "Dim"  "Neg"  "Neg"  "Dim"  "Dim"  "Dim"  "High" "Dim" 
##  [971] "Dim"  "High" "High" "Dim"  "Neg"  "Neg"  "Dim"  "High" "Dim"  "Dim" 
##  [981] "High" "Neg"  "High" "High" "Dim"  "Neg"  "Dim"  "High" "Neg"  "Neg" 
##  [991] "Neg"  "High" "High" "Dim"  "High" "Dim"  "High" "Dim"  "Neg"  "Dim" 
## [1001] "High" "Neg"  "Neg"  "Dim"  "Dim"  "High" "Neg"  "High" "High" "Dim" 
## [1011] "Neg"  "High" "Dim"  "Neg"  "Dim"  "Dim"  "Dim"  "Neg"  "Dim"  "Dim" 
## [1021] "Dim"  "Dim"  "Dim"  "High" "Dim"  "Dim"  "High" "Neg"  "High" "Neg" 
## [1031] "Neg"  "Neg"  "Dim"  "High" "High" "Neg"  "High" "Dim"  "Dim"  "Neg" 
## [1041] "Dim"  "Dim"  "Dim"  "High" "Dim"  "High" "Dim"  "Neg"  "Dim"  "Dim" 
## [1051] "High" "Dim"  "High" "Dim"  "Neg"  "High" "High" "Dim"  "Dim"  "Neg" 
## [1061] "Neg"  "Neg"  "Dim"  "High" "Dim"  "High" "Dim"  "Dim"  "Dim"  "High"
## [1071] "Dim"  "Neg"  "High" "Neg"  "Dim"  "Dim"  "Dim"  "High" "Dim"  "Neg" 
## [1081] "Neg"  "Neg"  "High" "Neg"  "High" "High" "High" "Dim"  "Dim"  "Neg" 
## [1091] "Dim"  "Neg"  "High" "Neg"  "Dim"  "Dim"  "Dim"  "Neg"  "Neg"  "Neg" 
## [1101] "Neg"  "High" "Dim"  "Dim"  "Dim"  "High" "Neg"  "Dim"  "Neg"  "Dim" 
## [1111] "Dim"  "High" "Neg"  "Dim"  "High" "Neg"  "High" "Neg"  "Dim"  "Dim" 
## [1121] "Dim"  "Neg"  "High" "Dim"  "Dim"  "Dim"  "High" "Dim"  "Dim"  "Dim" 
## [1131] "High" "Neg"  "Dim"  "Dim"  "Dim"  "Neg"  "Dim"  "Dim"  "Neg"  "High"
## [1141] "Dim"  "High" "Neg"  "Dim"  "Dim"  "Neg"  "High" "Dim"  "High" "Neg" 
## [1151] "Neg"  "Neg"  "High" "Neg"  "Dim"  "Dim"  "Dim"  "Dim"  "Dim"  "Neg" 
## [1161] "Neg"  "High" "High" "High" "High" "Dim"  "High" "Dim"  "Neg"  "Dim" 
## [1171] "Neg"  "Neg"  "Dim"  "High" "Dim"  "Neg"  "High" "Dim"  "High" "Dim" 
## [1181] "High" "High" "Dim"  "High" "High" "High" "Neg"  "Neg"  "High" "High"
## [1191] "Dim"  "Dim"  "Dim"  "Neg"  "Dim"  "Dim"  "Dim"  "Dim"  "Dim"  "Neg" 
## [1201] "High" "Neg"  "Dim"  "Dim"  "High" "High" "Dim"  "Neg"  "Dim"  "Dim" 
## [1211] "Dim"  "Dim"  "High" "High" "Neg"  "Dim"  "Dim"  "Dim"  "High" "High"
## [1221] "Dim"  "Dim"  "Dim"  "High" "High" "High" "Dim"  "Dim"  "Dim"  "High"
## [1231] "Neg"  "Dim"  "Neg"  "Dim"  "Dim"  "High" "High" "Dim"  "Dim"  "Dim" 
## [1241] "Neg"  "Dim"  "Dim"  "Dim"  "High" "Neg"  "Dim"  "Dim"  "Dim"  "Neg" 
## [1251] "High" "Dim"  "Dim"  "High" "High" "Dim"  "Dim"  "Neg"  "Dim"  "Neg" 
## [1261] "Neg"  "High" "Dim"  "Neg"  "Neg"  "Dim"  "Dim"  "Dim"  "High" "High"
## [1271] "Neg"  "High" "High" "Dim"  "High" "Neg"  "Dim"  "Dim"  "High" "Dim" 
## [1281] "High" "Dim"  "High" "Dim"  "Dim"  "Dim"  "High" "High" "Dim"  "Dim" 
## [1291] "Dim"  "Dim"  "Dim"  "High" "Dim"  "Neg"  "High" "Dim"  "Dim"  "Neg" 
## [1301] "Dim"  "Dim"  "Dim"  "High" "Dim"  "High" "Neg"  "High" "High" "Dim" 
## [1311] "Neg"  "Neg"  "Dim"  "Neg"  "Neg"  "Neg"  "Neg"  "Dim"  "Dim"  "Dim" 
## [1321] "Dim"  "High" "Dim"  "High" "Dim"  "High" "Dim"  "Dim"  "Dim"  "Dim" 
## [1331] "High" "Dim"  "Dim"  "Neg"  "Dim"  "Dim"  "High" "Dim"  "Dim"  "Dim" 
## [1341] "High" "Dim"  "Dim"  "Dim"  "Neg"  "Neg"  "Dim"  "High" "Neg"  "Dim" 
## [1351] "Neg"  "High" "High" "High" "Dim"  "Neg"  "Dim"  "Dim"  "Neg"  "Dim" 
## [1361] "Dim"  "Dim"  "Neg"  "Dim"  "Neg"  "Dim"  "Neg"  "Dim"  "High" "Dim" 
## [1371] "Dim"  "Dim"  "Neg"  "Dim"  "High" "High" "High" "Neg"  "Dim"  "Neg" 
## [1381] "Dim"  "High" "Neg"  "Dim"  "Dim"  "Neg"  "Neg"  "Neg"  "Dim"  "Dim" 
## [1391] "Dim"  "High" "Dim"  "Dim"  "Neg"  "Dim"  "Dim"  "Neg"  "Dim"  "Dim" 
## [1401] "High" "Dim"  "Neg"  "Dim"  "Dim"  "Neg"  "Dim"  "Dim"  "Neg"  "Dim" 
## [1411] "Dim"  "Neg"  "Dim"  "Dim"  "Dim"  "Neg"  "High" "Dim"  "High" "High"
## [1421] "Dim"  "Dim"  "Dim"  "Neg"  "Dim"  "High" "Neg"  "High" "Dim"  "High"
## [1431] "High" "Dim"  "Dim"  "Dim"  "Dim"  "Neg"  "High" "Neg"  "High" "Dim" 
## [1441] "Dim"  "Neg"  "Neg"  "Dim"  "Neg"  "Neg"  "Neg"  "High" "Neg"  "Dim" 
## [1451] "Dim"  "Dim"  "High" "Dim"  "Dim"  "High" "High" "Dim"  "Dim"  "Dim" 
## [1461] "Neg"  "High" "Neg"  "High" "Neg"  "High" "Neg"  "Neg"  "Dim"  "Dim" 
## [1471] "High" "Dim"  "Dim"  "Dim"  "Dim"  "Neg"  "Neg"  "Neg"  "Dim"  "Neg" 
## [1481] "Dim"  "High" "Dim"  "Dim"  "Neg"  "Neg"  "Neg"  "Dim"  "Neg"  "Dim" 
## [1491] "Dim"  "Dim"  "Dim"  "Dim"  "Neg"  "Neg"  "Dim"  "Dim"  "Neg"  "Neg" 
## [1501] "Dim"  "Neg"  "High" "Neg"  "High" "Dim"  "High" "High" "Neg"  "Dim" 
## [1511] "High" "Dim"  "Neg"  "Neg"  "Neg"  "Dim"  "Dim"  "Neg"  "High" "Dim" 
## [1521] "Dim"  "High" "High" "High" "Neg"  "Dim"

# Add ADT clust info to metadata
seurat.object.filt <- AddMetaData(object = seurat.object.filt, metadata = DNAM.clust, 
    col.name = "DNAM.clust.quantile")

head(seurat.object.filt@meta.data)

# Replace ident slot with DNAM clust info
clust.ident <- as.factor(seurat.object.filt@meta.data$DNAM.clust.quantile)
clust.ident

##    [1] Neg  Dim  High High Neg  Neg  Dim  Dim  Neg  Neg  Dim  Dim  Dim  High
##   [15] Dim  Dim  High Neg  Dim  Dim  Neg  Dim  Dim  Neg  Neg  High High Neg 
##   [29] Dim  High High Dim  Neg  High High Dim  High Neg  Dim  Neg  Dim  High
##   [43] Neg  Neg  Dim  High Neg  Dim  Dim  Dim  Dim  Dim  Neg  Dim  Dim  Dim 
##   [57] Dim  Neg  Dim  Dim  Dim  Dim  Dim  Dim  Dim  High Dim  Neg  Neg  Neg 
##   [71] Neg  Dim  Dim  High Dim  Dim  Neg  High High High High High High Dim 
##   [85] Dim  Dim  Neg  High High Dim  Dim  Dim  Dim  Dim  Dim  Dim  Neg  Dim 
##   [99] High High High Dim  Dim  High Dim  Dim  High Dim  Dim  Dim  Dim  High
##  [113] Dim  Neg  Neg  Dim  Dim  High Dim  Dim  High Dim  High High Dim  Dim 
##  [127] Dim  High Neg  High Neg  Dim  Dim  Dim  Dim  Dim  Neg  Dim  High Dim 
##  [141] Dim  Dim  High Neg  High Dim  Dim  Neg  Neg  Neg  Dim  Dim  High Neg 
##  [155] High Neg  Dim  Dim  Neg  Neg  High Dim  High Neg  Dim  Neg  Dim  Neg 
##  [169] Dim  Dim  Dim  Dim  Dim  Dim  Dim  Neg  High Dim  High Neg  High Neg 
##  [183] Dim  Dim  Neg  High High Neg  High Dim  Neg  Dim  Dim  High Neg  High
##  [197] Neg  Neg  Neg  Neg  High Dim  Dim  Neg  High Dim  Dim  Dim  Dim  High
##  [211] Neg  Neg  Dim  Dim  Dim  Dim  Dim  Neg  Dim  Dim  Neg  High Dim  Dim 
##  [225] Neg  Dim  High High Neg  High Dim  High Dim  High Dim  Neg  High Dim 
##  [239] High Dim  Neg  Neg  Dim  Dim  Dim  Dim  High Dim  Dim  Dim  High Dim 
##  [253] Dim  High Dim  Dim  Neg  Neg  Neg  Dim  High Neg  Dim  Dim  Dim  Dim 
##  [267] Neg  High Neg  High Dim  Neg  Neg  High Neg  Neg  Dim  High High Neg 
##  [281] Neg  Dim  High Neg  High Dim  High High Dim  Dim  Dim  High High Dim 
##  [295] Dim  High Dim  High High Dim  Neg  Dim  Dim  Dim  High Dim  Neg  Neg 
##  [309] Neg  High High Dim  Neg  Neg  Dim  Dim  Dim  High Dim  Dim  Dim  Dim 
##  [323] Neg  Neg  Dim  Dim  High Dim  Neg  Dim  High Neg  Dim  High High High
##  [337] High Dim  Dim  Neg  Dim  High Dim  Dim  High High Dim  Neg  Dim  Dim 
##  [351] Dim  Dim  Dim  Dim  Dim  Dim  Dim  Dim  High Neg  High Neg  Dim  Dim 
##  [365] Neg  Neg  Dim  Dim  Dim  Dim  Dim  Neg  Dim  High Dim  Dim  Dim  Dim 
##  [379] Dim  Neg  High Dim  Dim  Dim  High Dim  High Neg  Neg  Dim  Neg  Dim 
##  [393] Dim  High Dim  Neg  Neg  Dim  High Neg  Neg  Neg  Neg  Dim  Dim  Neg 
##  [407] Neg  High Dim  Dim  Dim  Dim  Neg  Neg  Dim  High Dim  High High High
##  [421] Dim  High Dim  Dim  Dim  Neg  Dim  High Dim  Neg  High Dim  Dim  Dim 
##  [435] Neg  Dim  Dim  Dim  Dim  High Dim  High High Neg  High Dim  High Neg 
##  [449] Dim  Dim  Neg  Dim  Neg  Dim  Dim  Dim  High Neg  High High Dim  Neg 
##  [463] High High Dim  Neg  Dim  Dim  Neg  Dim  Neg  Dim  Neg  High Dim  Neg 
##  [477] Dim  Dim  Dim  Dim  Dim  Dim  Dim  Dim  High High Neg  High Dim  Dim 
##  [491] Dim  Dim  Neg  Neg  High Dim  High Neg  Dim  Dim  Dim  Neg  High High
##  [505] Dim  Dim  High Dim  Neg  Neg  High Dim  Dim  High High Dim  Dim  Dim 
##  [519] Dim  Dim  High Neg  Dim  Neg  Neg  Dim  Neg  Dim  Dim  High Neg  Neg 
##  [533] Dim  Dim  Dim  Dim  High Dim  Neg  High Dim  Dim  High High Dim  Dim 
##  [547] High High Dim  Dim  Dim  High Dim  Dim  Dim  High Dim  High Dim  High
##  [561] Neg  Dim  Neg  Dim  Dim  Dim  Dim  High High Dim  Dim  Dim  High High
##  [575] Neg  High Neg  Dim  Neg  Dim  High High Dim  Neg  Dim  Dim  High High
##  [589] Neg  Neg  Dim  Dim  Dim  Neg  Dim  Dim  Dim  Dim  Neg  Neg  Neg  Neg 
##  [603] Dim  Dim  Dim  High Dim  Neg  Dim  Dim  High Neg  Dim  High Neg  Neg 
##  [617] Neg  Neg  Neg  High High High Dim  Neg  Neg  Dim  Dim  Dim  Dim  Dim 
##  [631] Dim  Dim  Neg  Dim  Dim  Dim  Dim  Dim  Dim  Neg  Neg  Neg  Dim  Dim 
##  [645] Neg  High High Dim  Neg  Neg  Dim  Neg  High Neg  Dim  Dim  Dim  High
##  [659] Dim  Neg  Dim  Dim  Dim  Dim  Dim  Dim  High Dim  High Dim  Neg  High
##  [673] Neg  Dim  Dim  Dim  Neg  High Neg  High Neg  Neg  Dim  Neg  Dim  Dim 
##  [687] Dim  Dim  Dim  High High High Dim  High Dim  High Dim  Neg  High Dim 
##  [701] Dim  Dim  Neg  Dim  Dim  High Neg  Dim  Dim  High Neg  High Neg  Dim 
##  [715] Dim  Dim  Dim  Dim  Dim  Dim  High Dim  Dim  Dim  High Dim  High Neg 
##  [729] Dim  High Dim  Neg  Neg  Dim  Dim  Neg  Dim  High Dim  Neg  Dim  Dim 
##  [743] Neg  Dim  High Dim  High Neg  Dim  High Dim  Dim  Dim  Neg  Dim  Dim 
##  [757] Neg  High Dim  Dim  Dim  Neg  Dim  High Dim  Dim  Neg  Neg  Dim  Neg 
##  [771] High Dim  Dim  Dim  High High Dim  Dim  Dim  Neg  High Dim  Neg  Neg 
##  [785] Dim  Dim  Dim  Neg  Dim  High Dim  High Dim  High High Dim  High Dim 
##  [799] Neg  Neg  Dim  Dim  Dim  Neg  Neg  Neg  Neg  High Dim  Neg  Dim  High
##  [813] Dim  Neg  Neg  Dim  Neg  Dim  Dim  Dim  High Neg  High Dim  Dim  Neg 
##  [827] High Neg  Neg  Dim  High Neg  High Dim  Neg  Dim  High High Neg  Dim 
##  [841] Dim  Dim  Dim  Dim  Dim  Dim  Dim  Neg  Dim  Neg  Dim  High Neg  Neg 
##  [855] Dim  Dim  High Neg  High High Dim  Dim  Dim  Neg  Neg  Neg  Dim  High
##  [869] High Dim  Dim  Dim  Dim  Dim  Dim  Dim  High High Dim  High Neg  Neg 
##  [883] High Neg  Dim  Neg  Dim  High High High Neg  High Dim  Dim  High Dim 
##  [897] Dim  Dim  Neg  High High Dim  Dim  Dim  Dim  Neg  High Neg  High Dim 
##  [911] Dim  High Dim  High Dim  Neg  Dim  Dim  High High High High High High
##  [925] High High High High Dim  Dim  Dim  Dim  Dim  Neg  Dim  Neg  Dim  High
##  [939] Dim  Dim  Dim  Neg  Neg  Dim  Dim  Dim  High Dim  High Neg  High Dim 
##  [953] High Neg  High Dim  Neg  Dim  Dim  Dim  Dim  High Dim  Neg  Neg  Dim 
##  [967] Dim  Dim  High Dim  Dim  High High Dim  Neg  Neg  Dim  High Dim  Dim 
##  [981] High Neg  High High Dim  Neg  Dim  High Neg  Neg  Neg  High High Dim 
##  [995] High Dim  High Dim  Neg  Dim  High Neg  Neg  Dim  Dim  High Neg  High
## [1009] High Dim  Neg  High Dim  Neg  Dim  Dim  Dim  Neg  Dim  Dim  Dim  Dim 
## [1023] Dim  High Dim  Dim  High Neg  High Neg  Neg  Neg  Dim  High High Neg 
## [1037] High Dim  Dim  Neg  Dim  Dim  Dim  High Dim  High Dim  Neg  Dim  Dim 
## [1051] High Dim  High Dim  Neg  High High Dim  Dim  Neg  Neg  Neg  Dim  High
## [1065] Dim  High Dim  Dim  Dim  High Dim  Neg  High Neg  Dim  Dim  Dim  High
## [1079] Dim  Neg  Neg  Neg  High Neg  High High High Dim  Dim  Neg  Dim  Neg 
## [1093] High Neg  Dim  Dim  Dim  Neg  Neg  Neg  Neg  High Dim  Dim  Dim  High
## [1107] Neg  Dim  Neg  Dim  Dim  High Neg  Dim  High Neg  High Neg  Dim  Dim 
## [1121] Dim  Neg  High Dim  Dim  Dim  High Dim  Dim  Dim  High Neg  Dim  Dim 
## [1135] Dim  Neg  Dim  Dim  Neg  High Dim  High Neg  Dim  Dim  Neg  High Dim 
## [1149] High Neg  Neg  Neg  High Neg  Dim  Dim  Dim  Dim  Dim  Neg  Neg  High
## [1163] High High High Dim  High Dim  Neg  Dim  Neg  Neg  Dim  High Dim  Neg 
## [1177] High Dim  High Dim  High High Dim  High High High Neg  Neg  High High
## [1191] Dim  Dim  Dim  Neg  Dim  Dim  Dim  Dim  Dim  Neg  High Neg  Dim  Dim 
## [1205] High High Dim  Neg  Dim  Dim  Dim  Dim  High High Neg  Dim  Dim  Dim 
## [1219] High High Dim  Dim  Dim  High High High Dim  Dim  Dim  High Neg  Dim 
## [1233] Neg  Dim  Dim  High High Dim  Dim  Dim  Neg  Dim  Dim  Dim  High Neg 
## [1247] Dim  Dim  Dim  Neg  High Dim  Dim  High High Dim  Dim  Neg  Dim  Neg 
## [1261] Neg  High Dim  Neg  Neg  Dim  Dim  Dim  High High Neg  High High Dim 
## [1275] High Neg  Dim  Dim  High Dim  High Dim  High Dim  Dim  Dim  High High
## [1289] Dim  Dim  Dim  Dim  Dim  High Dim  Neg  High Dim  Dim  Neg  Dim  Dim 
## [1303] Dim  High Dim  High Neg  High High Dim  Neg  Neg  Dim  Neg  Neg  Neg 
## [1317] Neg  Dim  Dim  Dim  Dim  High Dim  High Dim  High Dim  Dim  Dim  Dim 
## [1331] High Dim  Dim  Neg  Dim  Dim  High Dim  Dim  Dim  High Dim  Dim  Dim 
## [1345] Neg  Neg  Dim  High Neg  Dim  Neg  High High High Dim  Neg  Dim  Dim 
## [1359] Neg  Dim  Dim  Dim  Neg  Dim  Neg  Dim  Neg  Dim  High Dim  Dim  Dim 
## [1373] Neg  Dim  High High High Neg  Dim  Neg  Dim  High Neg  Dim  Dim  Neg 
## [1387] Neg  Neg  Dim  Dim  Dim  High Dim  Dim  Neg  Dim  Dim  Neg  Dim  Dim 
## [1401] High Dim  Neg  Dim  Dim  Neg  Dim  Dim  Neg  Dim  Dim  Neg  Dim  Dim 
## [1415] Dim  Neg  High Dim  High High Dim  Dim  Dim  Neg  Dim  High Neg  High
## [1429] Dim  High High Dim  Dim  Dim  Dim  Neg  High Neg  High Dim  Dim  Neg 
## [1443] Neg  Dim  Neg  Neg  Neg  High Neg  Dim  Dim  Dim  High Dim  Dim  High
## [1457] High Dim  Dim  Dim  Neg  High Neg  High Neg  High Neg  Neg  Dim  Dim 
## [1471] High Dim  Dim  Dim  Dim  Neg  Neg  Neg  Dim  Neg  Dim  High Dim  Dim 
## [1485] Neg  Neg  Neg  Dim  Neg  Dim  Dim  Dim  Dim  Dim  Neg  Neg  Dim  Dim 
## [1499] Neg  Neg  Dim  Neg  High Neg  High Dim  High High Neg  Dim  High Dim 
## [1513] Neg  Neg  Neg  Dim  Dim  Neg  High Dim  Dim  High High High Neg  Dim 
## Levels: Dim High Neg

seurat.object.filt <- SetIdent(seurat.object.filt, value = clust.ident)

head(seurat.object.filt@active.ident)

## AAACCCAGTAATCAAG AAACCCAGTCTGTCCT AAACGCTAGTCACGCC AAAGAACGTTAACAGA 
##              Neg              Dim             High             High 
## AAAGGATAGTCCTACA AAAGGATGTATCGTAC 
##              Neg              Neg 
## Levels: Dim High Neg

# Give Active ident and metadata columns factorisation with correct ordering
levels(seurat.object.filt@active.ident)

## [1] "Dim"  "High" "Neg"

seurat.object.filt@active.ident <- factor(seurat.object.filt@active.ident, levels = c("Neg", 
    "Dim", "High"))

levels(seurat.object.filt@active.ident)

## [1] "Neg"  "Dim"  "High"

levels(seurat.object.filt@meta.data$DNAM.clust.quantile)

## NULL

seurat.object.filt@meta.data$DNAM.clust.quantile <- factor(seurat.object.filt@meta.data$DNAM.clust.quantile, 
    levels = c("Neg", "Dim", "High"))

levels(seurat.object.filt@meta.data$DNAM.clust.quantile)

## [1] "Neg"  "Dim"  "High"

0.2.5 ADT cluster # and %

table(seurat.object.filt@meta.data$DNAM.clust)

## 
##  Neg  Dim High 
##  384  756  386

prop.table(table(seurat.object.filt@meta.data$DNAM.clust)) * 100

## 
##      Neg      Dim     High 
## 25.16383 49.54128 25.29489

0.2.6 Visualising ADT subsets

# Create output directory
if (!dir.exists("output/figures/ADT_subsetting")) {
    dir.create("output/figures/ADT_subsetting", recursive = T)
}

# Quick QC
VlnPlot(object = seurat.object.filt, features = c("nFeature_RNA", "nCount_RNA", "percent.mito"))

dev.copy(pdf, "output/QC/QC_nFeature_nCount_RNA_percentMito_Post_filtering_DNAM1_ADT_subsets.pdf")

## pdf 
##   3

dev.off()

## quartz_off_screen 
##                 2

# Overall ADT detection Using custom x axis

RidgePlot(seurat.object.filt, "adt_DNAM1", group.by = "orig.ident", log = TRUE, sort = FALSE, 
    cols = col.scheme) + ggplot2::coord_cartesian(xlim = c(0.9, 4))

dev.copy(pdf, "output/figures/ADT_subsetting/Ridgeplot_adt_DNAM1_total_log_cust_xaxis.pdf")

## pdf 
##   3

dev.off()

## quartz_off_screen 
##                 2

# Vis ADT clustering on VlnPlot
VlnPlot(seurat.object.filt, "adt_DNAM1", group.by = "DNAM.clust.quantile", sort = FALSE, 
    cols = col.scheme)

dev.copy(pdf, "output/figures/ADT_subsetting/Vlnplot_adt_DNAM1_groups_quantile.pdf")

## pdf 
##   3

dev.off()

## quartz_off_screen 
##                 2

# Vis ADT clustering on RidgePlot
RidgePlot(seurat.object.filt, "adt_DNAM1", group.by = "DNAM.clust.quantile", sort = FALSE, 
    cols = col.scheme)

dev.copy(pdf, "output/figures/ADT_subsetting/Ridgeplot_adt_DNAM1_groups_quantile.pdf")

## pdf 
##   3

dev.off()

## quartz_off_screen 
##                 2

# Visualise Cd226 gene expression

VlnPlot(seurat.object.filt, "Cd226", group.by = "DNAM.clust.quantile", sort = FALSE, 
    cols = col.scheme)

dev.copy(pdf, "output/figures/ADT_subsetting/Vlnplot_adt_DNAM1_groups_Cd226_expression_quantile.pdf")

## pdf 
##   3

dev.off()

## quartz_off_screen 
##                 2

RidgePlot(seurat.object.filt, "Cd226", group.by = "DNAM.clust.quantile", sort = FALSE, 
    cols = col.scheme)

dev.copy(pdf, "output/figures/ADT_subsetting/Ridgeplot_adt_DNAM1_groups_Cd226_expression_quantile.pdf")

## pdf 
##   3

dev.off()

## quartz_off_screen 
##                 2

0.3 Differential gene expression

0.3.1 Calculating DEGs

# DEG of High vs. Neg ADT subsets
DNAMhigh.vs.neg.markers.wilcox <- FindMarkers(seurat.object.filt, ident.1 = "High", 
    ident.2 = "Neg", logfc.threshold = 0.25, min.pct = 0.1, test.use = "wilcox", 
    only.pos = FALSE)

write.csv(DNAMhigh.vs.neg.markers.wilcox, "output/tables/DNAMhigh_vs_neg_markers_wilcox.csv")
print(paste0("Number of sig differentially expressed genes = ", sum(DNAMhigh.vs.neg.markers.wilcox$p_val_adj < 
    0.05)))

## [1] "Number of sig differentially expressed genes = 214"

# Do diff gene expression with no threshold or filtering
DNAMhigh.vs.Neg.markers.wilcox.NO.FILT <- FindMarkers(seurat.object.filt, ident.1 = "High", 
    ident.2 = "Neg", logfc.threshold = 0, min.pct = 0, test.use = "wilcox", only.pos = FALSE)

write.csv(DNAMhigh.vs.Neg.markers.wilcox.NO.FILT, "output/tables/DNAMhigh_vs_Neg_markers_wilcox_ALL_genes_NO_filt.csv")
print(paste0("Number of sig differentially expressed genes = ", sum(DNAMhigh.vs.Neg.markers.wilcox.NO.FILT$p_val_adj < 
    0.05)))

## [1] "Number of sig differentially expressed genes = 900"

0.3.2 Write data to file

save.data.frame.function <- function(df, title) {
    
    x <- as.data.frame(as.matrix(df))
    
    write.table(x, paste0("output/tables/", title, ".txt"), sep = "\t", quote = FALSE)
}

# RNA dataset Write raw data
save.data.frame.function(seurat.object.filt@assays$RNA@counts, "Raw_dataframe_RNA")

# write filtered data
save.data.frame.function(seurat.object.filt@assays$RNA@data, "filtered_dataframe_RNA")

# Write metadata
save.data.frame.function(seurat.object.filt@meta.data, "Meta_data_dataframe")

# ADT dataset Write raw data
save.data.frame.function(seurat.object.filt@assays$ADT@counts, "Raw_dataframe_ADT")

# write filtered data
save.data.frame.function(seurat.object.filt@assays$ADT@data, "filtered_dataframe_ADT")

0.4 Visualisation

0.4.1 Volcano DEGs

if(!dir.exists("output/figures/VolcanoPlots")){
  dir.create("output/figures/VolcanoPlots", 
             recursive = TRUE)}

#############
# Function
#############

clean.data <- function(x){
  
  # data = data.frame with cols (genes, PValue/FDR, logFC)
  
  colnames(x) <- c("PValue", "logFC", "pct.1", "pct.2", "FDR")
  colnames(x) 
  x$genes <- rownames(x)
  head(x)
  return(x)
  
}



# Format data
volcano.data <- clean.data(DNAMhigh.vs.Neg.markers.wilcox.NO.FILT)



# Create colour key for logFC and FDR sig value of interest
FDR.cutoff <- 0.05
logFC.cutoff <- 0.25

# Establish custom colour key for increased vs. decreased genes         
# set the base colour as 'black'
keyvals <- rep('grey50', nrow(volcano.data))

# set the base name/label as 'NS'
names(keyvals) <- rep('NS', nrow(volcano.data))

# modify keyvals for variables with fold change > 0.25 & FDR < 0.05
keyvals[which(volcano.data$logFC > logFC.cutoff & volcano.data$FDR < FDR.cutoff)] <- DNAM1.high.col
names(keyvals)[which(volcano.data$logFC > logFC.cutoff & volcano.data$FDR < FDR.cutoff)] <- 'High'

# modify keyvals for variables with fold change < -0.25 & FDR < 0.05
keyvals[which(volcano.data$logFC < -logFC.cutoff & volcano.data$FDR < FDR.cutoff)] <- DNAM1.neg.col
names(keyvals)[which(volcano.data$logFC < -logFC.cutoff & volcano.data$FDR < FDR.cutoff)] <- 'Neg'

unique(names(keyvals))

## [1] "High" "Neg"  "NS"

# Plot Volcano with top 10 DEGs and custom GOI

# top 10 up and down DEGs
sig.genes <- volcano.data[volcano.data$FDR < 0.05, ]

goi.top <- top_n(sig.genes, 10, logFC)
goi.bottom <- top_n(sig.genes, -10, logFC)

# Hand picked GOIs
curated.goi <- c("Cd226", "Eomes", "Rora", 
                 "Ifng", "Cd69", "Gzma", 
                 "Gzmb", "Cd48", "Gzmk", 
                 "Icam1", "Crtam", "Lag3", 
                 "Tnfrsf9", "Cd247")
curated.goi

##  [1] "Cd226"   "Eomes"   "Rora"    "Ifng"    "Cd69"    "Gzma"    "Gzmb"   
##  [8] "Cd48"    "Gzmk"    "Icam1"   "Crtam"   "Lag3"    "Tnfrsf9" "Cd247"

goi <- unique(c(goi.top$genes, goi.bottom$genes, curated.goi))

# Plot highly customised Volcano         
EnhancedVolcano(volcano.data,
                lab = rownames(volcano.data),
                x = 'logFC',
                y = 'FDR',
                selectLab = goi,
                xlim = c(-1.5, 2.6),
                ylim = c(0, 80),
                title = "DNAM1 High vs. Neg",
                subtitle = "",
                ylab = bquote(~-Log[10]~adjusted~italic(P)),
                pCutoff = 0.05,
                FCcutoff = 0.25,
                pointSize = 5, 
                #col = c("black", "black", "black", "red2"),
                colCustom = keyvals,
                colAlpha = 0.2, 
                drawConnectors = TRUE,
                widthConnectors = 0.5,
                colConnectors = 'grey30',
                typeConnectors = "open",
                labSize = 4,
                gridlines.major = TRUE,
                gridlines.minor = FALSE,
                border = 'partial',
                borderWidth = 0.5,
                borderColour = 'black')   

dev.copy(pdf, "output/figures/VolcanoPlots/Volcano_DNAM1_High_vs_Neg_top_10up_and_10down_Plus_curatedGOI.pdf")

## pdf 
##   3

dev.off()

## quartz_off_screen 
##                 2

# Please note 
# A number of genes were removed as part of figure generation. 
# Namely, AA467197, Gm12840, Slc17a6
# These were removed as they were deemed uninformative.

0.4.2 Heatmaps

# Create output directory
if (!dir.exists("output/figures/heatmaps")) {
    dir.create("output/figures/heatmaps", recursive = T)
}

# Subset dataset for easier plotting
seurat.object.filt.small <- subset(seurat.object.filt, downsample = 300)



# Top 20 Up- and down-regulated genes in High Vs. Neg comparison

sig.genes <- DNAMhigh.vs.neg.markers.wilcox[DNAMhigh.vs.neg.markers.wilcox$p_val_adj < 
    0.05, ]

sig.genes$GeneID <- rownames(sig.genes)

up.genes <- top_n(sig.genes, 20, avg_logFC)$GeneID
down.genes <- top_n(sig.genes, -20, avg_logFC)$GeneID

plot.genes <- c(up.genes, down.genes)

# Plot heatmap
DoHeatmap(seurat.object.filt.small, features = unique(plot.genes), assay = "RNA", 
    group.colors = col.scheme, angle = 90, raster = FALSE)

dev.copy(pdf, "output/figures/heatmaps/DEGs_top20UP_DN.pdf")

## pdf 
##   3

dev.off()

## quartz_off_screen 
##                 2

0.5 GSEA

0.5.1 Convert Symbol to EntrezID

# Using Bitr package
keytypes(org.Mm.eg.db)

##  [1] "ACCNUM"       "ALIAS"        "ENSEMBL"      "ENSEMBLPROT"  "ENSEMBLTRANS"
##  [6] "ENTREZID"     "ENZYME"       "EVIDENCE"     "EVIDENCEALL"  "GENENAME"    
## [11] "GO"           "GOALL"        "IPI"          "MGI"          "ONTOLOGY"    
## [16] "ONTOLOGYALL"  "PATH"         "PFAM"         "PMID"         "PROSITE"     
## [21] "REFSEQ"       "SYMBOL"       "UNIGENE"      "UNIPROT"

Entrez.ids <- bitr(rownames(DNAMhigh.vs.Neg.markers.wilcox.NO.FILT), fromType = "SYMBOL", 
    toType = "ENTREZID", OrgDb = "org.Mm.eg.db")

# 6.54% failed to map

# Merge dataset with pval and fold change info
GSEA.dataframe <- merge(Entrez.ids, DNAMhigh.vs.Neg.markers.wilcox.NO.FILT, by.x = "SYMBOL", 
    by.y = "row.names")

dim(GSEA.dataframe)  # 12,127 genes 

## [1] 12127     7

0.5.2 msigDB

# Analysis using msig database

# Msig Database gene sets: H: hallmark gene sets C1: positional gene sets C2:
# curated gene sets C3: motif gene sets C4: computational gene sets C5: GO gene
# sets C6: oncogenic signatures C7: immunologic signatures


# Isolate just the differentially expressed genes
sig.genes <- GSEA.dataframe[GSEA.dataframe$p_val_adj < 0.05, ]

dim(sig.genes)  # 869 genes

## [1] 869   7

# Subset by upregulated or downregulated genes
up.sig.genes <- sig.genes[sig.genes$avg_logFC > 0, ]
dim(up.sig.genes)  # 753 genes

## [1] 753   7

# Set variables
cat.val <- "C5"
pval.cut.val <- 0.05
qval.cut.val <- 0.2

# Formatting database
m_t2g <- msigdbr(species = "Mus musculus", category = paste0(cat.val)) %>% dplyr::select(gs_name, 
    entrez_gene)


# upregulated
e.up <- enricher(up.sig.genes$ENTREZID, TERM2GENE = m_t2g, pvalueCutoff = pval.cut.val, 
    qvalueCutoff = qval.cut.val)

e.up <- setReadable(e.up, org.Mm.eg.db, keyType = "ENTREZID")

print(head(e.up))

##                                                                                                            ID
## GO_T_CELL_ACTIVATION                                                                     GO_T_CELL_ACTIVATION
## GO_LEUKOCYTE_CELL_CELL_ADHESION                                               GO_LEUKOCYTE_CELL_CELL_ADHESION
## GO_LYMPHOCYTE_DIFFERENTIATION                                                   GO_LYMPHOCYTE_DIFFERENTIATION
## GO_T_CELL_DIFFERENTIATION                                                           GO_T_CELL_DIFFERENTIATION
## GO_REGULATION_OF_T_CELL_ACTIVATION                                         GO_REGULATION_OF_T_CELL_ACTIVATION
## GO_POSITIVE_REGULATION_OF_LEUKOCYTE_CELL_CELL_ADHESION GO_POSITIVE_REGULATION_OF_LEUKOCYTE_CELL_CELL_ADHESION
##                                                                                                   Description
## GO_T_CELL_ACTIVATION                                                                     GO_T_CELL_ACTIVATION
## GO_LEUKOCYTE_CELL_CELL_ADHESION                                               GO_LEUKOCYTE_CELL_CELL_ADHESION
## GO_LYMPHOCYTE_DIFFERENTIATION                                                   GO_LYMPHOCYTE_DIFFERENTIATION
## GO_T_CELL_DIFFERENTIATION                                                           GO_T_CELL_DIFFERENTIATION
## GO_REGULATION_OF_T_CELL_ACTIVATION                                         GO_REGULATION_OF_T_CELL_ACTIVATION
## GO_POSITIVE_REGULATION_OF_LEUKOCYTE_CELL_CELL_ADHESION GO_POSITIVE_REGULATION_OF_LEUKOCYTE_CELL_CELL_ADHESION
##                                                        GeneRatio   BgRatio
## GO_T_CELL_ACTIVATION                                      82/685 448/16881
## GO_LEUKOCYTE_CELL_CELL_ADHESION                           61/685 334/16881
## GO_LYMPHOCYTE_DIFFERENTIATION                             60/685 340/16881
## GO_T_CELL_DIFFERENTIATION                                 48/685 242/16881
## GO_REGULATION_OF_T_CELL_ACTIVATION                        52/685 315/16881
## GO_POSITIVE_REGULATION_OF_LEUKOCYTE_CELL_CELL_ADHESION    42/685 211/16881
##                                                              pvalue
## GO_T_CELL_ACTIVATION                                   1.854225e-31
## GO_LEUKOCYTE_CELL_CELL_ADHESION                        1.677063e-23
## GO_LYMPHOCYTE_DIFFERENTIATION                          2.466146e-22
## GO_T_CELL_DIFFERENTIATION                              2.865564e-20
## GO_REGULATION_OF_T_CELL_ACTIVATION                     3.564098e-18
## GO_POSITIVE_REGULATION_OF_LEUKOCYTE_CELL_CELL_ADHESION 6.108670e-18
##                                                            p.adjust
## GO_T_CELL_ACTIVATION                                   9.211790e-28
## GO_LEUKOCYTE_CELL_CELL_ADHESION                        4.165824e-20
## GO_LYMPHOCYTE_DIFFERENTIATION                          4.083938e-19
## GO_T_CELL_DIFFERENTIATION                              3.559030e-17
## GO_REGULATION_OF_T_CELL_ACTIVATION                     3.541288e-15
## GO_POSITIVE_REGULATION_OF_LEUKOCYTE_CELL_CELL_ADHESION 4.650126e-15
##                                                              qvalue
## GO_T_CELL_ACTIVATION                                   7.539865e-28
## GO_LEUKOCYTE_CELL_CELL_ADHESION                        3.409733e-20
## GO_LYMPHOCYTE_DIFFERENTIATION                          3.342709e-19
## GO_T_CELL_DIFFERENTIATION                              2.913072e-17
## GO_REGULATION_OF_T_CELL_ACTIVATION                     2.898550e-15
## GO_POSITIVE_REGULATION_OF_LEUKOCYTE_CELL_CELL_ADHESION 3.806136e-15
##                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                          geneID
## GO_T_CELL_ACTIVATION                                   Adam8/Anxa1/Apbb1ip/B2m/Bcl3/Camk4/Card11/Casp8/Ccl5/Ccnd3/Ccr2/Cd274/Cd44/Cd47/Cd80/Cd86/Cd8b1/Cdc42/Chd7/Coro1a/Ctla4/Dock2/Dock8/Egr1/Flot2/Gata3/Gpr18/Gpr183/Grap2/Gsn/H2-M3/Havcr2/Icam1/Icos/Ifng/Il18r1/Il2ra/Il7r/Irf1/Irf4/Itgal/Itk/Klrk1/Lck/Lcp1/Lfng/Lgals1/Lgals3/Lgals9/Lilr4b/Map3k8/March7/Msn/Nfkbid/Nfkbiz/Nlrc3/Pag1/Pik3cd/Ppp3ca/Ppp3cb/Prr7/Psmb10/Ptger4/Pycard/Rac2/Rasal3/Rasgrp1/Rhoa/Ripk2/Rora/Runx2/Runx3/Satb1/Sema4a/Smad7/Socs1/Spn/Thy1/Tnfaip8l2/Vsir/Zc3h12a/Zfp683
## GO_LEUKOCYTE_CELL_CELL_ADHESION                                                                                                                                                    Adam8/Anxa1/Card11/Ccl5/Ccr2/Cd274/Cd44/Cd47/Cd80/Cd86/Cdc42/Coro1a/Ctla4/Cx3cr1/Dock8/Ezr/Flot2/Gata3/Grap2/H2-M3/Havcr2/Icam1/Icos/Ifng/Il2ra/Il7r/Irf1/Itgal/Itgb1/Itgb2/Itgb7/Klrk1/Lck/Lgals1/Lgals3/Lgals9/Lilr4b/Map3k8/March7/Msn/Nfkbid/Nfkbiz/Pag1/Pycard/Rac2/Rasal3/Rasgrp1/Rhoa/Ripk2/Runx3/Selplg/Skap1/Smad7/Socs1/Spn/Stk10/Thy1/Tnfaip8l2/Vsir/Wnk1/Zc3h12a
## GO_LYMPHOCYTE_DIFFERENTIATION                                                                                                                                                             Adam8/Anxa1/B2m/Bcl3/Camk4/Card11/Ccr2/Cd80/Cd86/Chd7/Ctla4/Dock10/Dock2/Egr1/Flt3l/Gata3/Gpr18/Gpr183/H2-M3/Id2/Ifng/Il18r1/Il2ra/Il7r/Inpp5d/Irf1/Irf2bp2/Irf4/Itgb1/Itk/Klf6/Lck/Lfng/Lgals1/Lgals9/Lilr4b/Nfkbid/Nfkbiz/Notch2/Pglyrp2/Pik3cd/Plcl2/Ppp3cb/Prr7/Ptger4/Ptk2b/Rasgrp1/Rhoa/Ripk2/Rora/Runx2/Runx3/Satb1/Sema4a/Smad7/Socs1/Spn/Vsir/Zc3h12a/Zfp683
## GO_T_CELL_DIFFERENTIATION                                                                                                                                                                                                                                              Adam8/Anxa1/B2m/Bcl3/Camk4/Card11/Ccr2/Cd80/Cd86/Chd7/Ctla4/Dock2/Egr1/Gata3/Gpr18/Gpr183/H2-M3/Ifng/Il18r1/Il2ra/Il7r/Irf1/Irf4/Itk/Lck/Lfng/Lgals9/Lilr4b/Nfkbid/Nfkbiz/Pik3cd/Ppp3cb/Prr7/Ptger4/Rasgrp1/Rhoa/Ripk2/Rora/Runx2/Runx3/Satb1/Sema4a/Smad7/Socs1/Spn/Vsir/Zc3h12a/Zfp683
## GO_REGULATION_OF_T_CELL_ACTIVATION                                                                                                                                                                                                     Adam8/Anxa1/Camk4/Card11/Ccl5/Ccr2/Cd274/Cd47/Cd80/Cd86/Cdc42/Coro1a/Ctla4/Dock8/Flot2/Gata3/Grap2/Gsn/H2-M3/Havcr2/Icos/Ifng/Il2ra/Il7r/Irf1/Irf4/Klrk1/Lck/Lgals1/Lgals3/Lgals9/Lilr4b/Map3k8/March7/Nfkbid/Nfkbiz/Pag1/Pycard/Rac2/Rasal3/Rasgrp1/Rhoa/Ripk2/Runx3/Smad7/Socs1/Spn/Thy1/Tnfaip8l2/Vsir/Zc3h12a/Zfp683
## GO_POSITIVE_REGULATION_OF_LEUKOCYTE_CELL_CELL_ADHESION                                                                                                                                                                                                                                               Adam8/Anxa1/Card11/Ccl5/Ccr2/Cd274/Cd44/Cd47/Cd80/Cd86/Cdc42/Coro1a/Ctla4/Dock8/Flot2/Gata3/Grap2/H2-M3/Havcr2/Icam1/Icos/Ifng/Il2ra/Il7r/Klrk1/Lck/Lgals1/Lgals9/Lilr4b/Map3k8/Nfkbid/Nfkbiz/Pycard/Rasal3/Rasgrp1/Rhoa/Ripk2/Runx3/Skap1/Socs1/Thy1/Vsir
##                                                        Count
## GO_T_CELL_ACTIVATION                                      82
## GO_LEUKOCYTE_CELL_CELL_ADHESION                           61
## GO_LYMPHOCYTE_DIFFERENTIATION                             60
## GO_T_CELL_DIFFERENTIATION                                 48
## GO_REGULATION_OF_T_CELL_ACTIVATION                        52
## GO_POSITIVE_REGULATION_OF_LEUKOCYTE_CELL_CELL_ADHESION    42

# Visualisation of enrichment
print(barplot(e.up, showCategory = 10, font.size = 7) + ggtitle(paste0("Upregulated genes msigDB ", 
    cat.val[i], " enrichment")))

dev.copy(pdf, paste0("output/figures/GSEA/MsigDB_", cat.val[i], "_enrichment_barplot_Upregulated.pdf"))

## pdf 
##   3

dev.off()

## quartz_off_screen 
##                 2

# write files to disk
write.csv(e.up@result, paste0("output/tables/GSEA/MsigDB_", cat.val[i], "_enrich_upreg.csv"))

0.6 Imputation

0.6.1 Imputation calculation

# Impute gene expression
imputed.seurat <- RunALRA(seurat.object.filt, genes.use = rownames(seurat.object.filt))

# The matrix went from 20.75% nonzero to 42.98% nonzero

0.7 Imputed data visualisation

0.7.1 mRNA vs Protein scatter plots

if (!dir.exists("output/figures/scatterplot")) {
    dir.create("output/figures/scatterplot", recursive = T)
}


# Pearson correlation between the two features is displayed above the plot.

FeatureScatter(imputed.seurat, "Cd226", "adt_DNAM1", cols = col.scheme)

dev.copy(pdf, "output/figures/scatterplot/Cd226_vs_DNAM1_Imputed.pdf")

## pdf 
##   3

dev.off()

## quartz_off_screen 
##                 2

0.7.2 VlnPlots

if(!dir.exists("output/figures/VlnPlots_manuscript")){dir.create("output/figures/VlnPlots_manuscript", recursive = T)}


Genes.to.plot <- c("Cd69", "Cd226", "Crtam", "Cd27", "Cd28",
                   "Tnfrsf4", "Tnfrsf9", "Furin", "Il2ra", "Il12rb2",
                   "Il7r", "Ccr7", "Cd44", "Mki67", "Cd2",
                   "Ifng", "Nkg7", "Klrg1", "Slamf1", "Slamf7",
                   "Klrc1", "S1pr1", "Fasl", "Fas", "Cx3cr1",
                   "Ctsw", "Cd96", "Entpd1", "Pdcd1", "Adora2a",
                   "Lag3", "Tgfb1", "Havcr2", "Ctla4", "Cd160",
                   "Cd274", "Tigit", "Ikzf1", "Ikzf2", "Ikzf3",
                   "Tox", "Eomes", "Tbx21", "Id2", "Gata3",
                   "Irf4", "Irf8", "Rora", "Klf2", "Klf3", 
                   "Bach2", "Zeb2", "Myb", "Tcf7", "Itgb2",
                   "Cdc42", "Itgal", "Itgax", "Itgb1", "Cd48",
                   "Tln1", "Rap1a", "Rap1b", "Rassf5", "Rac2", 
                   "Ccl5", "Ccl4", "Xcl1", "Il16", "Csf1",
                   "Gzma", "Gzmb", "Gzmk")

for(i in 1:length(Genes.to.plot)){
  
  print(VlnPlot(imputed.seurat, 
                paste0(Genes.to.plot[i]),
                cols = col.scheme, 
                #ncol = 3,
                pt.size = 0))
  dev.copy(pdf, paste0("output/figures/VlnPlots_manuscript/VlnPlot_", Genes.to.plot[i], ".pdf"))
  dev.off()
  
}